function COut = StaVelMain(fileStations,OptAct,OptsGen)
%StaVelMain - General StaVel function
%
%COut = StaVelMain(fileStations,OptAct,OptsGen)
%
% This function, which calls the external package Hector, can performs 
% all computations necessary to obtain the velocities of the GNSS stations
% listed in the Excel or ASCII file fileStations. 
%
% fileStations
% ------------
% This file contains the names of all stations to be included in the 
% calculations.
% Each station must be represented by its standard 4-character name. 
% In the case of an Excel file, the station names must be placed in the 
% first column (no more than the first column is read). 
% In the case of an ASCII file, only a column is admitted. 
% If fileStations is undefined or empty, the filename can be interactively 
% managed.
%
% OptAct
% ------
% The second input argument, OptAct, allows the choice of the operation
% to be carried out:
%   1.  ALL STAGES, NO CME REMOVAL (TSDATA & RAW MOM FILE GENERATION, 
%           OFFSET SEARCH/TAKING, OUTLIER REMOVAL, TREND ESTIMATION)
%   2.  TSDATA AND RAW MOM FILE GENERATION ONLY
%   3.  OFFSET SEARCH/TAKING ONLY (TS/MOM FILES REQUIRED)
%   4.  OUTLIER REMOVAL ONLY (OFFSETS REQUIRED)
%   5.  TREND ESTIMATION ONLY (OUTLIERS REQUIRED)
%   6.  OUTLIER REMOVAL, TREND ESTIMATION (OFFSETS REQUIRED)
%   7.  CME ESTIMATION AND FILTERING (TRENDS REQUIRED)
%   8.  TREND RE-ESTIMATION FOR CME-FREE TIME SERIES (CME FILTERING REQUIRED)
%   9.  TRENDS ALREADY COMPUTED, OUTPUT CELL VARIABLE GENERATION ONLY
%   10. EXIT
% If OptAct is undefined, empty or outside 1:10, the choice can be carried 
% out in an interactive way.
%
% If OptAct is 1 or 3 (chosen as an input argument or by means of interactive
% choice), the option about the offset management can be interactively
% chosen:
%       1. automatic search based on Hector (not recommended);
%       2. manual search;
%       3. taken from Nevada Geodetic Laboratory database;
%       4. taken from NGL, then manually verified and integrated; 
%       5. taken from another database;
%       6. taken from another database, then manually verified and integrated.
%
% If OptAct is 7 (chosen as an input argument or by means of interactive
% choice), the option about the CME estimation method can be interactively
% chosen:
%   1:  simple stacking;
%   2:  stacking weighted by errors on daily positions;
%   3:  stacking weighted by distance of each station from the network 
%       center;
%   4:  stacking weighted by correlation between stations time series;
%   5:  simple stacking, with recalculation of possible existing 
%       detrended zero-mean time series;
%   6:  stacking weighted by errors on daily positions, with recalculation 
%       of possible existing detrended zero-mean time series;
%   7:  stacking weighted by distance of each station from the network 
%       center, with recalculation of possible existing detrended zero-mean
%       time series;
%   8:  stacking weighted by correlation between stations time series, with 
%       recalculation of possible existing detrended zero-mean time series;
%   9:  simple stacking with minimal detrending, i.e. detrending based on 
%       simple least square fitting to a straight line (in this case, the 
%       Hector-based trend is not used);
%   10.	stacking weighted by errors on daily positions with minimal 
%       detrending;
%   11.	stacking weighted by distance of each station from the network 
%       center with minimal detrending;
%   12.	stacking weighted by correlation between stations time series with 
%       minimal detrending.
%
% If OptAct is 9 (chosen as an input argument or by means of interactive
% choice) and CME-filtered data are available, the use can generate the
% output variable COut on the basis or non filtered or filtered data. The
% option can be interactively chosen:
%   1: original trends;
%   2: trends obtained from CME-filtered time series.
%
% If OptAct is a 2-length vector, OptAct(1) is the option 1-10 as above, 
% and OptAct(2) is the option for the offset management, provided that
% OptAct(1) is 1 or 3, the option for the CME estimation, provided that
% OptAct(1) is 7, or the option for the choice of original or CME-filtered
% trends, provided that OptAct(1) is 9 and these trends exist. The second 
% element, i.e. OptAct(2), is actually used only OptAct(1) is 1, 3, 7 or 9.
% The vector option on OptAct is introduced to allow the use of StaVelMain 
% under GNU Octave.
%
% OptsGen
% -------
% OptsGen is the struct variable with the general StaVel options, generated
% by means of geneOpts function. If OptsGen is undefined or empty, the name 
% of the file carrying such a struct variable can be interactively managed.
% In this case, the file must have a field whose name is either Opts or
% OptsGen.
%
% Output
% ------
% The output variable COut is non-empty if OptAct is 1, 5, 6, 8 or 9, empty 
% elsewere. If COut is non-empty, it is a N-by-11 cell array, where N is
% the number of managed stations, and the columns are:
% 
%   [station name, initial time, final time, station latitude, 
%       station longitude, station height (meters a.s.l.), 
%       estimated trend East, estimated trend East standard deviation,
%       estimated trend North, estimated trend North standard deviation,
%       estimated trend Up, estimated trend East standard deviation Up]
%
% If OptAct is 8 or [9 2], the trends in COut are CME-filtered.
%
% See also tsData, geneOpts, CMEstackFiltering, InspectOffsetSearch. 

% G. Teza, 2022.

if nargin < 3
    OptsGen = [];
end
if nargin < 2
    OptAct = [];
end
if nargin < 1
    fileStations = [];
end

lastOptAct = 10;

numOptAct = length(OptAct);     % for the use of Octave 
if numOptAct == 2
    IOptAct = false;             % to avoid the option about OptOff
    if ismember(OptAct(1),[1 3]) 
        OptOff = OptAct(2);
    elseif OptAct(1) == 7
        OptStack = OptAct(2);
    elseif OptAct(1) == 9
        OptTrend = OptAct(2);
    end
    OptAct = OptAct(1);
else
    OptStack = [];
    OptTrend = [];
    IOptAct = true;
end

if isempty(OptAct) || ~ismember(OptAct,1:10)
    figIn = figure('Name','StaVel main function');
    set(figIn,'Units','normalized','OuterPosition',[0 0 1 1]);
    imshow([])
    title('StaVel general option','FontSize',24,'Color','b');  
    axis equal;
    lsgen = {...
        '1.  ALL STAGES, NO CME REMOVAL (TSDATA/RAW MOM FILE GENERATION, OFFSETS SEARCH/TAKING, OUTLIERS REMOVAL, TREND ESTIMATION)',...
        '2.  TSDATA AND RAW MOM FILE GENERATION ONLY',...
        '3.  OFFSETS SEARCH/TAKING ONLY (TSDATA/RAW MOM FILES REQUIRED)',...
        '4.  OUTLIERS REMOVAL ONLY (OFFSETS REQUIRED)',...
        '5.  TREND ESTIMATION ONLY (OUTLIERS REQUIRED)',...
        '6.  OUTLIER REMOVAL, TREND ESTIMATION (OFFSETS REQUIRED)',...
        '7.  CME ESTIMATION AND FILTERING (TRENDS REQUIRED)',...
        '8.  TREND RE-ESTIMATION FOR CME-FILTERED TIME SERIES (CME FILTERING REQUIRED)',...
        '9.  TRENDS ALREADY COMPUTED, OUTPUT CELL VARIABLE GENERATION ONLY',...
        '10. EXIT'};
    lbIn = uicontrol(figIn,'Style','listbox',...
        'string',lsgen,...
        'FontSize',14,...
        'Units','normalized','Position', [0.08 0.3 0.88 0.36],...
        'Callback', 'uiresume(gcbf)');
    uiwait(figIn);
    OptAct = lbIn.Value;
    delete(lbIn);
    close(figIn);
    if OptAct == 10
        OptAct = lastOptAct;
    end
end
if ismember(OptAct,[1 3]) && IOptAct
    figOff = figure('Name','Offset option');
    set(figOff,'Units','normalized','OuterPosition',[0 0 1 1]);
    imshow([])
    title('Offsets search option','FontSize',24,'Color','b');  
    axis equal;
    lsOff = {...
        '1. AUTOMATIC OFFSETS SEARCH',...
        '2. MANUAL OFFSETS SEARCH',...
        '3. OFFSETS FROM NEVADA GEODETIC LABORATORY (NGL) DATABASE',...
        '4. OFFSETS FROM NGL, THEN MANUALLY CHECKED/INTEGRATED',...
        '5. OFFSETS FROM ANOTHER DATABASE',...
        '6. OFFSETS FROM ANOTHER DATABASE, THEN MANUALLY CHECKED/INTEGRATED',...
        '7. EXIT'};
     lbOff = uicontrol(figOff,'Style','listbox',...
        'String',lsOff,...
        'FontSize',14,...
        'Units','normalized','Position', [0.1 0.3 0.8 0.40],...
        'Callback', 'uiresume(gcbf)');
     uiwait(figOff);
     OptOff = lbOff.Value;
     delete(lbOff);
     close(figOff);
     if OptOff == 7  % exit command
        OptAct = lastOptAct;
     end
end
if ismember(OptAct,9) && isempty(OptTrend)
    figT = figure('Name','Trend option');
    set(figT,'Units','normalized','OuterPosition',[0 0 1 1]);
    imshow([])
    title('Trend data option','FontSize',24,'Color','b');  
    axis equal;
    lsT = {...
        '1. GENERATE OUTPUT CELL ARRAY WITH TRENDS FROM ORIGINAL DATA',...
        '2. GENERATE OUTPUT CELL ARRAY WITH TRENDS FROM CME-FILTERED DATA',...
        '3. EXIT'};
    lbT = uicontrol(figT,'Style','listbox',...
        'String',lsT,...
        'FontSize',14,...
        'Units','normalized','Position', [0.1 0.3 0.8 0.30],...
        'Callback', 'uiresume(gcbf)');
    uiwait(figT);
    OptTrend = lbT.Value;
    delete(lbT);
    close(figT);
    if OptTrend == 3  % exit command
        OptAct = lastOptAct;
    end
end
if ismember(OptAct,[1 3])  
    if ismember(OptOff,[3 4])
        UrlOrFilena = [];
        OffsetsString = 'Offsets taken from Nevada Geodetic Laboratory database';
    elseif ismember(OptOff,[5 6])
        [filena,pathna,findex] = uigetfile(...
            '*.mat', 'PLEASE SELECT THE OFFSETS DATABASE FILE');
        if findex
            UrlOrFilena = fullfile(pathna,filena);
            OffsetsString = sprintf('Offsets taken from file %s',UrlOrFilena); 
        else
            OptAct = lastOptAct;
            fprintf('Invalid choice, program stop \n');
        end
    end
end
if OptAct == lastOptAct
    COut = [];
    return
end

[b,fileStations] = readStationList(fileStations);
if isempty(b) 
    COut = [];
    return
else
    nin = size(b,1);
end

OptsGen = OptsGenLoad(OptsGen);
if isempty(OptsGen)
    COut = [];
    return
end

% the first row is assumed to be the header
COut = cell(nin,12);

try     % (for compatibility with a previous version)
    dirEnv = OptsGen.dirEnv;
catch
    dirEnv = OptsGen.dirNeu;
end
dirEnv = testDir(dirEnv);
dirTs  = testDir(OptsGen.dirTs);
dirRaw = testDir(OptsGen.dirRaw);
dirObs = testDir(OptsGen.dirObs);
dirPre = testDir(OptsGen.dirPre);
dirMom = testDir(OptsGen.dirMom);
dirPreCME = testDir(OptsGen.dirPreCME);
dirMomCME = testDir(OptsGen.dirMomCME);
dirCtl = testDir(OptsGen.dirCtl);
dirOut = testDir(OptsGen.dirOut);
dirJSON = testDir(OptsGen.dirJSON);
extTenv = OptsGen.extTenv;
dirHec = OptsGen.dirHector;
comps  = OptsGen.components;
try     % optional additional string for Env files
    AddEnv = OptsGen.AddEnv;
catch
    AddEnv = OptsGen.AddNeu;
end
AddTs  = OptsGen.AddTs;     % optional additional string for Ts files 
AddRaw = OptsGen.AddRaw;    % optional additional string for Raw files
AddObs = OptsGen.AddObs;    % optional additional string for Obs files
AddPre = OptsGen.AddPre;    % optional additional string for Pre files
AddMom = OptsGen.AddMom;    % optional additional string for Mom files
AddPreCME = OptsGen.AddPreCME;    % optional additional string for Pre CME-free files
AddMomCME = OptsGen.AddMomCME;    % optional additional string for Mom CME-free files
% AddOut = OptsGen.AddOut;    % optional additional string for Out files
AddJSON = OptsGen.AddJSON;  % optional additional string for JSON files
Lmin = OptsGen.Lmin;
if isempty(Lmin)
    Lmin = 0;
end

OptsOffset = OptsGen.Offsets;
OptsOffset.DataDirectory = OptsGen.dirRaw;
fileInitOffset = fullfile(dirCtl,'findoffset.ctl');

OptsOutlier = OptsGen.Outliers;
OptsOutlier.DataDirectory = OptsGen.dirObs;
fileInitOutlier = fullfile(dirCtl,'removeoutliers.ctl');

OptsTrend = OptsGen.Trend;
OptsTrend.DataDirectory = OptsGen.dirPre;
fileInitTrend = fullfile(dirCtl,'estimatetrend.ctl');

compsc = {'e','n','v'};    % please note the order "Hector style" 

OffsetsCell = [];   % to have non more than a database load
CME = [];           % to have no more than a CME computation  
for k = 1:nin
    namek = b{k,1};
    namek = strtrim(namek);     % to remove possible leading and trailing spaces
    namek = upper(namek);       % to warrant station name with capital letters
    if isempty(namek)
        continue
    end
    
    % tsData/mom file generation:
    ITs = false;
    ILok = 0;
    if ismember(OptAct,1:2)
        filenak = fullfile(dirEnv,[namek AddEnv extTenv]);
        if exist(filenak,'file')
            ts = tsData(filenak);
            Is = searchFilledFields(ts);
            if sum(Is) > 0  % at least a field is non-empty: neu/mom files can be generated
                if Lmin > 0
                    t = ts.t;
                    ti = t(1);
                    tf = t(end-1);
                    if (tf-ti)/365.25 > Lmin
                        ITs = true;
                        ILok = 1;
                    else
                        ILok = 2;
                    end
                else
                    ITs = true;
                    ILok = 1;
                end
                if ITs
                    save(fullfile(dirTs,[namek AddTs '.mat']),'ts');
                    fprintf('station %s - tsData file generated \n',namek);
                end
            end
        end
        if ~ITs
            if ILok == 0 
                fprintf('station %s - invalid data - no tsData file generated \n',namek);
            else
                fprintf('station %s - too short time series - no tsData file generated \n',namek);
            end
            continue
        end
    
        % raw mom file generation:
        if ITs 
            nameOutMomk = fullfile(dirRaw,[namek AddRaw]);
            namets = ts;                                    % tsData object
            [ITMom,~] = writeNeuMom(nameOutMomk,1,namets);
            fprintf('station %s - Raw mom files generated \n',namek);
        else
            ITMom = false;
        end
    end
    
    if (OptAct == 2) || ((OptAct == 1) && ~ITMom)
        continue
    end

    fileTs = fullfile(dirTs,[namek AddTs '.mat']); % defined once and for all

    vComp = OptsGen.components;
    if numel(vComp) == 3
        hin  = 1;
        hfin = 3;
    elseif numel(vComp) == 2
        hin  = 1;
        hfin = 2;
    else
        hin  = 3;
        hfin = 3;
    end
    
    for h = hin:hfin
        if ismember(h-1,comps)
            comph = compsc{h};
            
            % OFFSETS RECOGNITION OR OFFSETS ACQUISITION FROM DATABASE ----
            if ismember(OptAct,[1 3])
                stcOffkh = 1;
                
                if OptOff == 1
                    
                    ts = tsDataFileIn(fileTs);
                    if isempty(ts)
                        continue
                    end
                    % Automatic offset recognition
                    nameInkh = [namek AddRaw '_' num2str(h-1)];
                    nameOutkh = [namek AddObs '_' num2str(h-1)];
                    nameJSONkh = [namek AddJSON '_' num2str(h-1)];
                    fileJSONOffsetkh = fullfile(dirJSON,['findoffset_' nameJSONkh '.json']);
                    OptsOffsetkh = OptsOffset;
                    OptsOffsetkh.DataFile = [nameInkh '.mom'];
                    if exist(fullfile(dirRaw,OptsOffsetkh.DataFile),'file') ~= 2
                        continue
                    end
                    OptsOffsetkh.OutputFile = [dirObs '/' nameOutkh '.mom']; % here, fullfile cannot be used!!
                    fileFinOffsetkh  = [dirCtl '/findoffset_' nameOutkh '.ctl'];
                    stfOffkh = ctlFileEdit(fileInitOffset,fileFinOffsetkh,OptsOffsetkh);
                    if stfOffkh
                        fprintf('Station %s, component %d (%s) - offset recognition in progress... \n',namek,h-1,upper(comph));
                        strCommandOffsetkh = ['wsl ' dirHec 'findoffset ' fileFinOffsetkh];
                        [stcOffkh,resOffkh] = system(strCommandOffsetkh);
                        if stcOffkh == 0     % successful offset computation
                            txtFileOffkh = fullfile(dirOut,['Res_Offset_' nameOutkh '.txt']);
                            fidOffkh = fopen(txtFileOffkh,'w+');
                            fprintf(fidOffkh,resOffkh);
                            fclose(fidOffkh);
                            % JSON file management and JSON -> tsData transfer
                            movefile('findoffset.json',fileJSONOffsetkh);
                            if h == 1
                                offsetsE = jsonRead(fileJSONOffsetkh);
                                ts.offsetsE = offsetsE;
                            elseif h == 2
                                offsetsN = jsonRead(fileJSONOffsetkh);
                                ts.offsetsN = offsetsN;
                            else
                                offsetsV = jsonRead(fileJSONOffsetkh);
                                ts.offsetsV = offsetsV;
                            end
                            if h == hin
                                ts.offsetsOpts = OptsGen.Offsets;
                            end
                            save(fileTs,'ts')                            % to solve a possible Hector's problem
                            [headhk,datahk] = txtReadData(OptsOffsetkh.OutputFile);
                            if size(datahk,2) == 3     % in this case, the third column must be removed
                                datahk = datahk(:,1:2);
                                txtWriteData(OptsOffsetkh.OutputFile,headhk,datahk);
                            end
                            fprintf('Station %s, component %d (%s) - offset recognized \n',namek,h-1,upper(comph));
                        else
                            fprintf('Station %s, component %d (%s) - unsuccessful offset recognition \n',namek,h-1,upper(comph));
                        end
                    end

                elseif OptOff == 2

                    % Manual offset recognition
                    if h == hin
                        [~,ts] = InspOffsetMom(namek,OptsGen,vComp);
                        ts.offsetsOpts = 'Manual offsets recognition';
                        save(fileTs,'ts');
                    end
                    stcOffkh = 0;
    
                else

                    % Offsets from a database
                    if isempty(OffsetsCell)
                        fprintf('Offsets database loading in progress...\n');
                        OffsetsCell = NevadaAllOffsets(UrlOrFilena);
                        fprintf('... Offsets database loaded \n');
                    end
                    if h == hin
                        [IDOK,~,ts] = OffsetsfromCell(namek,OptsGen,OffsetsCell);
                        if ~isempty(ts)
                            if IDOK
                                fprintf('Station %s, offsets from database taken \n',namek);
                            else
                                fprintf('Station %s, no offsets from database \n',namek);
                            end
                            ts.offsetsOpts = OffsetsString;
                            save(fileTs,'ts');
                            if ismember(OptOff,[4 6])
                                [~,ts] = InspOffsetMom(namek,OptsGen,vComp,1);
                                ts.offsetsOpts = 'Data from Database manually checked/integrated';
                                save(fileTs,'ts');
                            end
                        else
                            fprintf('Station %s, no ts file available \n',namek);
                        end
                    end
                    stcOffkh = 0;

                end
                
                if (OptAct == 3) || (ismember(OptAct,[1 2]) && stcOffkh > 0)
                    continue
                end

            end
            
            % OUTLIER REMOVAL ---------------------------------------------
            if ismember(OptAct,[1 4 6])

                stcOutlierkh = 1;
                ts = tsDataFileIn(fileTs);
                if isempty(ts)
                    continue
                end
                nameInkh = [namek AddObs '_' num2str(h-1)];
                nameOutkh = [namek AddPre '_' num2str(h-1)];
                nameJSONkh = [namek AddJSON '_' num2str(h-1)];
                fileJSONOutlierkh = fullfile(dirJSON,['removeoutliers_' nameJSONkh '.json']);
                OptsOutlierkh = OptsOutlier;
                OptsOutlierkh.DataFile = [nameInkh '.mom'];
                if exist(fullfile(dirObs,OptsOutlierkh.DataFile),'file') ~= 2
                    continue
                end
                OptsOutlierkh.OutputFile = [dirPre '/' nameOutkh '.mom']; % here, fullfile cannot be used!!
                fileFinOutlierkh  = [dirCtl '/removeoutliers_' nameOutkh '.ctl'];
                stfOutlierkh = ctlFileEdit(fileInitOutlier,fileFinOutlierkh,OptsOutlierkh);
                if stfOutlierkh
                    fprintf('station %s, component %d (%s) - outlier removal in progress...\n',namek,h-1,upper(comph));
                    strCommandOutlierskh = ['wsl ' dirHec 'removeoutliers ' fileFinOutlierkh];
                    [stcOutlierkh,resOutlierkh] = system(strCommandOutlierskh);
                    if stcOutlierkh == 0     % successful outlier removal computation
                        txtFileOutlierkh = fullfile(dirOut,['Res_Outlier_' nameOutkh '.txt']);
                        fidOutlierkh = fopen(txtFileOutlierkh,'w+');
                        fprintf(fidOutlierkh,resOutlierkh);
                        fclose(fidOutlierkh);
                        % JSON file management and JSON -> tsData transfer
                        movefile('removeoutliers.json',fileJSONOutlierkh);
                        % to acquire for ts the cleaned data
                        [~,datahk] = txtReadData(OptsOutlierkh.OutputFile);
                        thk = MJD2serial(datahk(:,1));
                        phk = datahk(:,2);  % pre-processed (in mm!)
                        if h == 1
                            outliersE = jsonRead(fileJSONOutlierkh);
                            ts.outliersE = outliersE;
                            ts.cleanedE.t = thk;
                            ts.cleanedE.p = phk/1000; % to have data in m
                        elseif h == 2
                            outliersN = jsonRead(fileJSONOutlierkh);
                            ts.outliersN = outliersN;
                            ts.cleanedN.t = thk;
                            ts.cleanedN.p = phk/1000; % to have data in m
                        else
                            outliersV = jsonRead(fileJSONOutlierkh);
                            ts.outliersV = outliersV;
                            ts.cleanedV.t = thk;
                            ts.cleanedV.p = phk/1000; % to have data in m
                        end
                        ts.outliersOpts = OptsOutlier;
                        save(fileTs,'ts');
                        fprintf('station %s, component %d (%s) - outlier removed \n',namek,h-1,upper(comph));
                    else
                        fprintf('station %s, component %d (%s) - unsuccessful outlier removal \n',namek,h-1,upper(comph));
                    end
                end
                if (OptAct == 4) || (stcOutlierkh > 0)
                    continue
                end
            end

            % CME ESTIMATION AND FILTERING (OPTIONAL) ---------------------
            
            if OptAct == 7
                if isempty(CME)
                    if (hin == 1) && (hfin == 3)
                        NC = 3;
                    elseif (hin == 1) && (hfin == 2)
                        NC = 2;
                    else
                        NC = [];
                    end
                    CME = CMEstackFiltering(fileStations,OptsGen,OptStack,NC);
                end
            end
            
            % TREND COMPUTATION -------------------------------------------
            if ismember(OptAct,[1 5 6 8])
                
                stcTrendkh = 1;
                ts = tsDataFileIn(fileTs);
                if isempty(ts)
                    continue
                end
                if OptAct < 8
                    dirPree = dirPre;
                    AddPree = AddPre;
                    dirMome = dirMom;
                    AddMome = AddMom;
                    Adde = '';
                else
                    dirPree = dirPreCME;
                    AddPree = AddPreCME;
                    dirMome = dirMomCME;
                    AddMome = AddMomCME;
                    Adde = 'CME';
                end
                nameInkh = [namek AddPree '_' num2str(h-1)];
                nameOutkh = [namek AddMome '_' num2str(h-1)];
                nameJSONkh = [namek Adde AddJSON '_' num2str(h-1)];
                fileJSONTrendkh = fullfile(dirJSON,['estimatetrend_' nameJSONkh '.json']);
                OptsTrendkh = OptsTrend;
                OptsTrendkh.DataDirectory = dirPree;
                OptsTrendkh.DataFile = [nameInkh '.mom'];
                if exist(fullfile(dirPree,OptsTrendkh.DataFile),'file') ~= 2
                    continue
                end
                OptsTrendkh.OutputFile = [dirMome '/' nameOutkh '.mom']; % here, fullfile cannot be used!!
                fileFinTrendkh  = [dirCtl '/estimatetrend_' nameOutkh '.ctl'];
                stfTrendkh = ctlFileEdit(fileInitTrend,fileFinTrendkh,OptsTrendkh);
                if stfTrendkh
                    fprintf('station %s, component %d (%s) - trend estimation in progress... \n',namek,h-1,upper(comph));
                    [stcTrendkh,resTrendkh] = system(['wsl ' dirHec 'estimatetrend ' fileFinTrendkh]);
                    if stcTrendkh == 0     % successful trend computation
                        txtFileTrendkh = fullfile(dirOut,['Res_Trend_' nameOutkh '.txt']);
                        fidTrendkh = fopen(txtFileTrendkh,'w+');
                        fprintf(fidTrendkh,resTrendkh);
                        fclose(fidTrendkh);
                        % JSON file management and JSON -> tsData transfer
                        ts = tsDataFileIn(fileTs);
                        movefile('estimatetrend.json',fileJSONTrendkh);
                        [~,datahk] = txtReadData(OptsTrendkh.OutputFile);
                        thk = MJD2serial(datahk(:,1));
                        ohk = datahk(:,2);  % observation (in mm!)
                        mhk = datahk(:,3);  % model (in mm!)
                        if h == 1
                            estimatedTrendE = jsonRead(fileJSONTrendkh);
                            if OptAct < 8
                                ts.estimatedTrendE.trend = estimatedTrendE;
                                ts.estimatedTrendE.t = thk;
                                ts.estimatedTrendE.o = ohk/1000;
                                ts.estimatedTrendE.m = mhk/1000;
                            else
                                ts.CMEfiltered.estimatedTrendE.trend = estimatedTrendE;
                                ts.CMEfiltered.estimatedTrendE.t = thk;
                                ts.CMEfiltered.estimatedTrendE.o = ohk/1000;
                                ts.CMEfiltered.estimatedTrendE.m = mhk/1000;
                            end
                        elseif h == 2
                            estimatedTrendN = jsonRead(fileJSONTrendkh);
                            if OptAct < 8
                                ts.estimatedTrendN.trend = estimatedTrendN;
                                ts.estimatedTrendN.t = thk;
                                ts.estimatedTrendN.o = ohk/1000; % to have data in m
                                ts.estimatedTrendN.m = mhk/1000;
                            else
                                ts.CMEfiltered.estimatedTrendN.trend = estimatedTrendN;
                                ts.CMEfiltered.estimatedTrendN.t = thk;
                                ts.CMEfiltered.estimatedTrendN.o = ohk/1000;
                                ts.CMEfiltered.estimatedTrendN.m = mhk/1000;
                            end
                        else
                            estimatedTrendV = jsonRead(fileJSONTrendkh);
                            if OptAct < 8
                                ts.estimatedTrendV.trend = estimatedTrendV;
                                ts.estimatedTrendV.t = thk;
                                ts.estimatedTrendV.o = ohk/1000;
                                ts.estimatedTrendV.m = mhk/1000;
                            else
                                ts.CMEfiltered.estimatedTrendV.trend = estimatedTrendV;
                                ts.CMEfiltered.estimatedTrendV.t = thk;
                                ts.CMEfiltered.estimatedTrendV.o = ohk/1000;
                                ts.CMEfiltered.estimatedTrendV.m = mhk/1000;
                            end
                        end
                        ts.estimatedTrendOpts = OptsTrend;
                        save(fileTs,'ts');
                        fprintf('station %s, component %d (%s) - trend estimated \n',namek,h-1,upper(comph));
                        fprintf('station %s - %s, %s tsData file updated \n',namek,upper(comph),fileTs); 
                    else
                        fprintf('station %s, component %d (%s) - unsuccessful trend estimation \n',namek,h-1,upper(comph));    
                    end
                end
                if stcTrendkh > 0
                    continue
                end
            end

            % OUTPUT CELL VARIABLE GENERATION -----------------------------
            if ismember(OptAct,[1 6 8 9])
                if (OptAct == 8) || ((OptAct == 9) && isequal(OptTrend,2))
                    OptD = false;   % if OptD is false, CME-filtered data are used
                else
                    OptD = true;    % if OptD is true, original data are used
                end
                ts = tsDataFileIn(fileTs);
                if isempty(ts)
                    continue
                end
                if OptD && isempty(ts.estimatedTrendE)
                    if h == hin
                        fprintf('station %s, no available trend data \n',namek);
                    end
                    continue
                end
                if ~OptD && (isempty(ts.CMEfiltered) || length(fieldnames(ts.CMEfiltered)) < 8) 
                    if h == hin
                        fprintf('station %s, no available trend from CME-filtered data \n',namek);
                    end
                    continue
                end
                if h == hin
                    COut(k,1) = {ts.statName};
                    t = ts.t;
                    COut(k,2) = {t(1)};
                    COut(k,3) = {t(end)};
                    COut(k,4) = {ts.statLat};
                    COut(k,5) = {ts.statLon};
                    COut(k,6) = {ts.statHeight};
                end
                if h == 1
                    if OptD
                        COut(k,7) = {ts.estimatedTrendE.trend.trend};
                        COut(k,10) = {ts.estimatedTrendE.trend.trend_sigma};
                    else
                        COut(k,7) = {ts.CMEfiltered.estimatedTrendE.trend.trend};
                        COut(k,10) = {ts.CMEfiltered.estimatedTrendE.trend.trend_sigma};
                    end
                elseif h == 2
                    if OptD
                        COut(k,8) = {ts.estimatedTrendN.trend.trend};
                        COut(k,11) = {ts.estimatedTrendN.trend.trend_sigma};
                    else
                        COut(k,8) = {ts.CMEfiltered.estimatedTrendN.trend.trend};
                        COut(k,11) = {ts.CMEfiltered.estimatedTrendN.trend.trend_sigma};
                    end
                elseif h == 3
                    if OptD
                        COut(k,9) = {ts.estimatedTrendV.trend.trend};
                        COut(k,12) = {ts.estimatedTrendV.trend.trend_sigma};
                    else
                        COut(k,9) = {ts.CMEfiltered.estimatedTrendV.trend.trend};
                        COut(k,12) = {ts.CMEfiltered.estimatedTrendV.trend.trend_sigma};
                    end
                end
                if h == hfin
                    fprintf('station %s, output cell row generated \n',namek);
                end
            end
                    
        end
    end
end

if ~ismember(OptAct,[1,5,6,8,9])
    COut = [];
end

%% ANCILLARY FUNCTIONS

function dirOut = testDir(dirIn)
% This function tests if dirIn, if non-empty, exists as a directory.
% If dirIn is non-empty and does not exist, the new directory dirIn is
% made.
% If dirIn is non-empty and exists, dirOut = dirIn.
% If dirIn is empty, dirOut = dirIn.
if ~isempty(dirIn) && (exist(dirIn,'dir') ~= 7)
    mkdir(dirIn);
end
dirOut = dirIn;